$onText
DISE-2024 (Dynamic Integrated Space-Economy) model
Scenario: No debris (Clean Space Environment)
Growth model with satellites
File: DISE_2024_00.gms
Bongers, A. and Torres, J. L. (2025). Optimal path for orbital debris (DISE-2024)

$offText

*   Simulation period
Set  t  'Time periods'                          /2023*2272/;

Sets
    tfirst(t)   'First period (2023)'
        tlast(t)    'Last period (T)'
        tnotlast(t) 'All periods but last' ;

    tfirst(t)$(ord(t)=1)       =   yes;
    tlast(t)$(ord(t)=card(t))  =   yes;
    tnotlast(t)                =   not tlast(t);

Parameters

*   Economic parameters

**  Preferences parameters
    
    rho         'Social intertemporal preference rate'   /0.015/
    sigma       'Relative risk-aversion'                /1.5/
    beta0       'Discount factor'
    gN0         'labor growth rate parameter'          /0.05/
    Nstar       'Asymptotic population'                /10200/
    gA0         'TFP growth rate'                      /0.015/
    deltaA      'TFP growth rate trend'                /0.001/
    gQ0         'Satellite ISTC growth rate'           /0.030/
    deltaQ      'Satellite ISTC grwoth rate trend'     /0.005/

**  Launch cost
    b0      'Initial proportion of launch cost'     /0.3/
    gb0     'Initial change in launch cost'         /0.0/
    deltab  'Launch cost growth rate trend'         /0.0/

**  Initial values economic variables

   	k0	    'Initial Earth capital'			/555.6987/         		
    	s0     	'Initial Satellite capital'			/1.1959/  			
    	              
    i0      'Initial capital investment'    
    h0  'Initial satellite investment'
    y0  'Initial output'                /184.65/
    ii0     'Initial total investment'
    c0  'Initial consumption'
    l0  'Initial value of launches'
    x0      'Initial value of satellites destroyed'
    N0      'Initial population'                /8056/
    
** Preference parameters


**  Technological parameters
      
    alpha1  'Output-capital elasticity'                 /0.3479/
    alpha2  'Output-satellite elasticity'               /0.0021/ 
    deltak  'Capital depreciation rate'                 /0.0700/
    deltas  'Satellite depreciation rate'               /0.1500/ 
    a0  'Initial value of TFP'                    
    q0      'Initial value of satellite ISTC'           /1.0/

*   Physical parameters

    eta 'Number of satellites per launch'       /13.6/
    mu  'Mapping parameter'
    
**  Initial values physical variables
    
    ss0 'Operational satellites 2023'           /8500/
    ;
*   Calculus

    mu  =   ss0/s0;
    a0  =   y0/(k0**alpha1*s0**alpha2*(N0)**(1-alpha1-alpha2));
 

Parameters
        beta(t) 'Discount factor'
        tval(t) 'numerical value of t';
    tval(t) =   ord(t)-1;



Parameters
        b(t)    'Fraction of launch cost'
        gb(t)   'Decline in the launch cost';
        
        gb(t)           =   gb0*exp(-deltab*(tval(t)-1));
        b("2023")       =   b0;
        loop(t,b(t+1)   =   b(t)*exp(gb(t)););

*	Technology

Parameters
*   Technology
        A(t)    'TFP'
        Q(t)    'Satellite ISTC'
        gA(t)   'Growth rate of TFG'
        gQ(t)   'Growth rate of Satellite ISTC'
        ;
    
        gA(t)       =   gA0*exp(-deltaA*(tval(t)-1));
        gQ(t)       =   gQ0*exp(-deltaQ*(tval(t)-1));
        a("2023")       =   a0;
        q("2023")       =   q0;
        loop(t,A(t+1)   =   a(t)*exp(gA(t)););
        loop(t,Q(t+1)   =   q(t)*exp(gQ(t)););
    
*   Population
Parameters
    N(t)    'Population';
    N("2023")       =   N0;
    loop(t,N(t+1)   =   N(t)*(Nstar/N(t))**gN0;);
	
* The terminal weight beta(tlast) computation
    	beta(tnotlast(t))	=	power(1+rho,-tval(t));
	beta(tlast(t)) 		= 	(1/rho)*power(1+rho,1-tval(t));
	    
Variables
	SU      'Total summatory discounted utility'
	y(t)	'Output'
   	c(t)    'Consumption'
   	ratcy(t)'Ratio consumption-output'
   	ik(t)	'Invesment in Earth capital'
   	is(t)	'Investment in space capital'
   	h(t)    'Launch cost'
   	k(t)    'Earth capital stock'
	s(t)	'Stock of satellites'
	l(t)	'Launches'
	ii(t)   'Investment'
	hh(t)	'Number of new satellites inserted into orbit'
	U(t)	'Instantaneous utility'
	DU(t)   'Discounted utility'
	ss(t)   'Number of operational satellites (stock)'
	ll(t)	'Number of launches'
	srk(t) 'Investment rate in Earth capital'
    srs(t)  'Investment rate in Space capital'
	;

Equations

* 	Economic variables

    	utility(t)      'Instantaneous utility'
    	sumdutility	    'Summatory of discounted utility'
    	production(t)   'Cobb-Douglas production function'
    	allocation(t)   'Household choose between consumption, investment and abatement'
	totalinvest(t)	'Total investment'
    	accumulak(t) 	'Capital accumulation'
	accumulas(t)	'Satellite accumulation'
    	ifinal(t)       'Minimal capital investment in final period'
	hfinal(t)	'Minimal satellite investment in final period'  
	iifinal(t)	'Minimal total investment in final period'
	ratiocy(t)      'Consumption-output ratio'
	ratioiy(t) 'Investment rate in Earth capital'
    ratiohy(t)  'Investment rate in Space capital'
    launchcost(t)   'Launch cost'

*	Physical variables

	launches(t)	'Launches'
	satellites(t)	'Number of operational satellites'
	newsat(t)	'Number of new satellites inserted into orbit'
	nlaunches(t)   'Number of launches'
	;


* 	Equations of the model

**	Utility

	utility(t)..	U(t)	=e=	((c(t)/N(t))**(1-sigma)-1)/(1-sigma);
	sumdutility..   SU  	=e=	sum(t, beta(t)*U(t)*N(t));

**	Production, consumption and investment

	production(t)..    		y(t) 		=e=    a(t)*(k(t)**alpha1)*(s(t)**alpha2)*((N(t))**(1-alpha1-alpha2));
	allocation(t)..    		c(t) 		=e=    y(t)-ik(t)-is(t);
	totalinvest(t)..		ii(t)		=e=	   ik(t)+is(t);
	accumulak(tnotlast(t)).. 	k(t+1) 		=l= (1-deltak)*k(t)+ik(t);
	accumulas(tnotlast(t))..	s(t+1)		=l=	(1-deltas)*s(t)+q(t)*(1-b(t))*is(t);
*ifinal(tlast)..             sum(t$tl(t+1), inv(t+1)/inv(t)-y("c",t+1)/y("c",t) =e= 0;
	ifinal(tlast)..     ik(tlast) 	=g= 	(0.02+deltak)*k(tlast);
	hfinal(tlast)..		is(tlast)	=g=	(0.02+deltas)*s(tlast)/q(tlast);
	iifinal(tlast)..	ii(tlast)	=e=	ik(tlast)+is(tlast);
	ratiocy(t)..		ratcy(t)    =e= c(t)/y(t);
    ratioiy(t)..            srk(t)      =e= ik(t)/y(t);
    ratiohy(t)..            srs(t)      =e= is(t)/y(t);
    launchcost(t)..         h(t)    =e= (1-b(t))*is(t);

**	Satellite sector

	launches(t)..	   l(t)	    =e=	is(t)-h(t);
	nlaunches(t)..	   ll(t)	=e=	mu*h(t)/eta;
	newsat(t)..	       hh(t)	=e=	eta*ll(t);
	satellites(t)..	   ss(t)	=e=	mu*s(t)/q(t);

* Bounds.
k.lo(t) = 0.001;
s.lo(t)	= 0.001;
c.lo(t) = 0.001;
y.lo(t) = 0.001;

* Initial conditions
k.fx(tfirst)    =   k0;
s.fx(tfirst)    =   s0;
y.fx(tfirst)    =   y0;


model DISE / all /;

solve DISE maximizing SU using nlp;
solve DISE maximizing SU using nlp;

* Export to Excel using GDX utilities
execute_unload "soldise00.gdx"
execute '=gdx2xls soldise00.gdx';